library(smashr)

Now we consider estimating a spatially structured mean \(\mu\) = (\(\mu_1\), . . . , \(\mu_2\) ) from Poisson data:

\[y_t ∼ Pois(\mu_t), (t = 1, . . . , T)\] # Simulation 1

n = 1000
t = 1:n
spike.f = function(x) (5 * exp(-100 * round((x - 100)/10)^2) + 5 * exp(-100 * round((x - 400)/10)^2) + 5 * exp(-100 * round((x - 800)/10)^2) + 2 * exp(-100 * round((x - 500)/100)^2) 
   )
mu.s = spike.f(t)

# POISSON CASE
# ------------
# Scale the signal to be non-zero and to have a low average intensity.
mu.t = 1 + mu.s

# Simulate an example dataset.
set.seed(101)
X.s = rpois(n, mu.t)

# Run smash.
mu.est = smash(X.s, "poiss")

# Plot the estimated mean function (red) against the ground-truth (black).
plot(mu.t, type = "l", col="darkgreen", lwd=1.5, ylab= "# mutations/window", xlab = "Genomic position (kb)", ylim=c(0,max(X.s)))
points(X.s,col="darkgreen", cex=0.2)
lines(mu.est, col = "salmon", lwd=1.5)

Here, green is the real parameter and data. red is the esitmated mean by smash.

Simulation 2

n = 500
t = 1:n
spike.f = function(x) ( 2 * exp(-100 * round((x - 100)/10)^2) + 2 * exp(-100 * round((x - 450)/10)^2) + 0.5 * exp(-100 * round((x - 300)/100)^2) 
   )
mu.s = spike.f(t)

# POISSON CASE
# ------------
# Scale the signal to be non-zero and to have a low average intensity.
mu.t = 1 + mu.s

# Simulate an example dataset.
for (i in 1:50){
set.seed(i)
X.s = rpois(n, mu.t)

# Run smash.
mu.est = smash(X.s, "poiss")

# Plot the estimated mean function (red) against the ground-truth (black).
print(i)
plot(mu.t, type = "l", col="darkgreen", lwd=1.5, ylab= "# mutations/window", xlab = "Genomic position (kb)", ylim=c(0,max(X.s)))
points(X.s,col="darkgreen", cex=0.2)
lines(mu.est, col = "salmon", lwd=1.5)}
## [1] 1
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 2

## [1] 3

## [1] 4

## [1] 5

## [1] 6

## [1] 7

## [1] 8

## [1] 9

## [1] 10

## [1] 11

## [1] 12
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 13
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 14

## [1] 15
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 16
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 17
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 18

## [1] 19

## [1] 20
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 21

## [1] 22
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 23

## [1] 24

## [1] 25

## [1] 26

## [1] 27

## [1] 28

## [1] 29

## [1] 30
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 31

## [1] 32

## [1] 33

## [1] 34

## [1] 35

## [1] 36

## [1] 37
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 38

## [1] 39

## [1] 40

## [1] 41
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## [1] 42

## [1] 43

## [1] 44

## [1] 45

## [1] 46

## [1] 47

## [1] 48

## [1] 49

## [1] 50

Here, green is the real parameter and data. red is the esitmated mean by smash.

Apply to real data

We use 20kb window size around region KCNQ2 (chr20:62,037,542).

region <- head(read.table("table.100k.window.stand.diff.cutoff.4.topwindow.across.all.mutation.type.CSV", sep=",", skip=1, stringsAsFactors = F), n=5)
dat <- read.table("table.ASCWGS_20180504.WGS1902_hg19_controls_SNV_remove_recurrent_mutation.20kwindow.csv")
for (i in 1:nrow(region)){
  dat0 <- dat[dat$chr==region[i,2] & dat$start > region[i,3] - 6e6 & dat$start < region[i,4] + 6e6, ]
  # Run smash.
  mu.est = smash(dat0$obs, "poiss")
  plot(dat0$start/1e6, dat0$obs, type = "l", col=NA, lwd=1.5, ylab= "# mutations/window", xlab = "Genomic position (Mb)", ylim=c(0,max(dat0$obs)), main = paste0(region[i,2],":",region[i,3]))
  points(dat0$start/1e6, dat0$obs, col="darkgreen", cex=0.2)
  lines(dat0$start/1e6, mu.est, col = "salmon", lwd=1.5)
}
## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.

## Warning in estimate_mixprop(data, g, prior, optmethod = optmethod, control
## = control, : Optimization failed to converge. Results may be unreliable.
## Try increasing maxiter and rerunning.